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Abstract 

The gradual insertion method for direct calculation of the chemical potential by molecular 
simulation is applied in the NpT ensemble to different quadrupolar two-centre Lennard- 
Jones fluids at high density state points. The results agree well with Widom's test particle 
insertion but show at very high densities signiflcantly smaller statistical uncertainties. 
The gradual insertion method, which is coupled here with preferential sampling, extends 
the density range where reliable information on the chemical potential can be obtained. 
Application details are reported. 
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1 Introduction 



All molecular simulation techniques for the calculation of the chemical potential which are 
accompanied with trial insertion (or deletion) of real or ghost particles become inaccurate 
and finally fail in the case of high densities, especially when applied to molecular fluids with 
strong interactions. This is the case for the conventional test particle insertion method by 
Widom [T] as well as for the Gibbs ensemble (GE) [2] and grand canonical Monte-Carlo 
ensemble (GG) [3], which demand insertion and deletion attempts of real particles. 

Therefore, different Monte-Garlo techniques have been proposed to improve the efficiency 
of the calculation of the chemical potential. Recent developments are, e. g., the expanded 
ensemble methods the particle deletion scheme [H], the augmented grand canonical 
ensemble [7j, or the scaled particle Monte-Garlo method [8j. An overview is given by Kofke 
and Gummings [9]. 

In the present work, an expanded ensemble method is used. A first version was published 
by Vorontsov-Velyaminov et al. [3]. That method was applied for gradual insertion by 
Nezbeda and Kolafa |5j , who additionally introduced preferential sampling. Related work 
with expanded ensemble methods is described e.g. in p^ll|12] . 




We use the gradual insertion method as proposed by Nezbeda and Kolafa [5] which is 
based on the following idea: instead of inserting or deleting a real particle, a fluctuating 
particle is introduced, that undergoes changes in a set of different discrete 'sizes' or states 
of coupling with all other particles. Additionally, the preferential sampling of the fluid in 
the vicinity of the fluctuating particle is introduced for a further increase of accuracy. 

This method was applied to different ensembles (NVT j5j, GG GE and NpT 
[T^fT3] ) using different interaction potentials such as the hard sphere fluid [S], the square- 
well fluid [T^, the spherical Lennard- Jones fluid [EUS!, and to ternary mixtures of hard 




2 



spheres and hard heteronuclear diatomics [T7] . 

In this work we apply the gradual insertion method for the first time to a polar in- 
termolecular potential, the two-centre Lennard- Jones plus point quadrupole fluid. That 
potential allows to describe the thermodynamic properties of many real fluids, like carbon 
dioxide, ethyne, propyne, or carbon disulfide with good accuracy, as has been shown in 
a comprehensive study |18|19] . The present study is performed in the NpT ensemble at 
high densities. The results are critically compared to test particle insertion. Details on 
the gradual insertion method, which is used here together with preferential sampling, are 
given. 

2 Method 

In the following we briefly summarize the gradual insertion method in a NVT ensemble, 
for more details see [5]. The key element of the method is to introduce a 'fluctuating 
particle' which can appear at different states of coupling with all other particles. This set 
of states starts from a pointwise, completely interactionless or decoupled state, up to the 
fully coupled state, in which the fluctuating particle has the same properties as the other 
particles. Each of the partially coupled states is related to a NVT sub-ensemble with no 
physical meaning. These sub-ensembles can be depicted in a scheme as follows: 

[N] ^ [N + ^,] ^ ... ^ [N + ^k-i] ^ [N + ^k], 

Wo Wi Wk-l Wk 

where [N+ipj] denotes a sub-ensemble with regular particles and one fluctuating particle 
in the state j, with j = 0, . . . , k. The interaction energy of the fluctuating particle with 
the remaining particles is denoted by ipj, in the fully coupled state by ipk, and in the fully 
decoupled state it is ipo = 0. Note that [N + ip^.] equals [A^ + 1]. The arrows describe the 
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transitions between neighbouring states. Finally, Wj denotes the weight of the state j. 

Additionally to the standard Monte-Carlo moves in a NVT ensemble, trial translation and 
trial rotation, the trial change of fluctuating particle state is necessary. The probability 
of accepting a change of the fluctuating particle from state i to state j is given by 



where 11 (j — >i) and 11 (i — >j) denote a priori transition probabilities of the respective 
changes and (3 is the Boltzmann factor. 

The chemical potential in the NVT ensemble is obtained from [S] 



where Prob [A^] and Prob [A^ + 1] are the probabilities to observe an ensemble with A^ and 
A^ + 1 particles, respectively. The paranthesis <> denote the ensemble average. Further- 
more, /i*'^(T) is the temperature dependent ideal part of the chemical potential. If the 
fluctuating particle is in the fully decoupled state, an insertion attempt into a random 
new position follows. When the fluctuating particle reaches the fully coupled state another 
particle is chosen randomly, to be treated as the fluctuating one. 

Improvement of efficiency often can be obtained by using the preferential sampling method, 
see e. g. Allen and Tildesley [3]. The idea is to sample particles preferentially in the vicinity 
of the fluctuating particle. Therefore, a function /(r) is introduced, denoting the proba- 
bility of an attempt to move a host particle in a distance r from the fluctuating particle. 
The following steps are necessary to keep the microreversibility condition [5]: 

(i) random choice of a host particle, 

(ii) acceptance of this choice with the probability f{roid), 





(3) 
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(iii) random choice of a new position and orientation of the host particle, 

(iv) acceptance of this choice with the probabihty min(l, f{rnew)/fijoid))i 

(v) the new configuration is accepted with the probabihty min(l, exp{— /?([/„e^ — Uoid)}), 

where U denotes the configurational energy. 

In order to extend the gradual insertion method to the NpT ensemble the volume fluctu- 
ation step has to be considered additionally. To set the pressure p, a volume change from 
Void to Vnew IS acccpted with the probability [3] 



This step was applied at random position within the Markov chain, independent from the 
state of the fluctuating particle. Similar to the NVT ensemble, only configurations with 
a fully coupled fluctuating particle contribute to the Monte-Carlo sampling of physical 
properties. 

As the volume fluctuates in the NpT ensemble, Eq. (3) has to be modified to claculate 
the chemical potential 



Up to our knowledge, no formal derivation of Eq. (5) has been given in the literature so 
far. It is not within the scope of this work to close that gap. We will limit ourselves to 
show that the results from Eq. (5) are in agreement with Widom's method, cf. section 4. 

3 Investigated model and technical details 

In the present investigation we consider the two-centre Lennard- Jones plus pointquadru- 
pole fluid (2CLJQ). It is composed of two identical Lennard- Jones sites a distance L 
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Pacc{Vold ^Vnew) = mlu (l, {Vnew /Vold)'^ " exp {-l3{p{Vnew - Void) + Unew - Uold)}) ■ (4) 





apart plus an ideal point quadrupole of moment Q placed in the geometric centre of the 
molecule |20j. The charges of the quadrupole are arranged along the molecular axis in the 
symmetric sequence +,—,—, + or, having the same energetic effect in pure fluids, — , +, 
+, — . A detailed description of this fluid is given in [2T] . 

The Lennard- Jones parameters a and e of the 2CLJQ pair potential were used for the 
reduction of the thermodynamic properties and the model parameters L and Q: temper- 
ature T* = kT/e, pressure p* = pcr^/e, density p* = pa^, configurational energy u* = u/e, 
elongation L* = L/cr, squared quadrupolar moment Q*"^ = Q"^/ (ecx^). 

For the simulation runs = 512 particles were used, the cut-off radius was set to 4 a, 
applying periodic boundary conditions and the minimum image convention. The long 
range corrections for the two-centre Lennard- Jones potential [22] were considered, in the 
case of quadrupolar interactions no long range corrections have to be made. The number 
of loops, defined below, was 50 000. The maximum values of translation distance, rotation 
angle, and volume change were adjusted to yield acceptance rates of roughly 0.5. Statistical 
uncertainties were calculated by conventional block averaging [23] • 

The number of states of the fluctuating particle was chosen to A; = 6. For the fluctuating 
particle in the state j we set cTj/a = (1/2 + ^J j/24:), ej/e = j/k, Lj/L = 1, and Q'^/Q^ = 
j/k. The diameter of the hard sphere, shielding the quadrupolar interaction site, was 
set to 0.4 a, except for the fully decoupled state j = 0, where it is zero. The shielding 
does not allow quadrupolar interaction sites to approach closer to each other than 0.4 a. 
Finally, the weights Wj had to be chosen. Starting from unity for all states, they were 
adjusted during an equilibration period in order to achieve a roughly equal distribution 
in the occurrence of all states. 

One Monte-Carlo loop is defined here as N trial translations, 2/3- N trial rotations, and 1 
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trial volume change, which are the regular NpT moves. The additional gradual insertion 
moves within a loop are: Mq-N attempts to change the state of the fluctuating particle, 
Mp-N attempts to translate the fluctuating particle, 2/3-Mf-N attempts to rotate the 
fluctuating particle, Mp-N attempts to find a host particle for preferential translation, 
and 2/3-Mp-N attempts to find a host particle for preferential rotation. Here, Mc, Mp, 
and Mp are the parameters of the gradual insertion method. 

4 Results and discussion 

In Fig. 1 running averages of the chemical potential for the 2CLJQ fluid (L* = 0.8, 
Q*2 = 4) in the liquid phase at T/T^ ^ 0.55 are shown function of Monte-Carlo loops, 
or respectively, molecular dynamics time steps. The technical details of the molecular 
dynamics simulations are given in [21]. The Monte-Carlo loops, as defined here, and the 
molecular dynamics time steps should be roughly comparable in the sense that both result 
in a new configuration of the entire system. The better convergence of gradual insertion at 
this high density state point can clearly be seen. The large steps in the Widom curves are 
typical when this method is applied to very dense fluids. They are due to the bad statistics, 
which comes from the small number of test particles which contribute significantly to the 
average. 

In Table 1, a comparison is given between simulation data from the present work and 
results from Widom's method [T] taken from ^T\. Six different 2CLJQ fluids are investi- 
gated in the liquid state close to the bubble line at two temperatures (T/T^ ~ 0.55 and 
T/Tc ~ 0.8, respectively). The first finding is, that the results of both methods agree 
within their statistical uncertainties, which is a numerical proof of Eq. (5). Regarding, 
secondly, the statistical uncertainties for all 2CLJQ fluids at the higher temperatures 
(T/Tc ~ 0.8), where intermediate liquid densities are found, both methods yield results 
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with approximately the same statistics. But for T/Tp f« 0.55, where the densities are 
roughly 20% higher, the uncertainties of Widoms's method increase by a factor of 10 to 
15. On the other hand, the uncertainties of the gradual insertion increase only by a fac- 
tor of about 2 to up to 4. So, at state points with high densities (T/Tc ^ 0.55) gradual 
insertion yields results that have considerably better statistics. For the 2CLJQ fluid with 
the highest elongation and highest quadrupole (L* = 0.8, Q*^ = 4), at the low tempera- 
ture, Widom's method is close to failure (the uncertainty in absolute numbers is almost 
±1). Gradual insertion gives a decent result with ±0.02. In such very dense and srongly 
interacting fluids, the accuracy of Widom's method can not be improved significantly by 
using more test particles or by increasing the length of the simulation run. 

In order to compare the computational effort, we performed Monte-Carlo simulations in 
the NpT ensemble on a conventional workstation applying both methods and keeping all 
other parameters, such as number of particles or cut-off radius constant. It turns out, that 
the simulation run with gradual insertion {Mc=10, Mp=10, Mp=50) needs 19.7 CPU h 
and the run with Widom's method (2 000 test particles after each loop) needs 3.7 CPU h. 
Therein, regular NpT configuration generation consumes 1.6 CPU h. So the computational 
effort for gradual insertion is an order of magnitude higher. For a comparison with the 
same CPU time spent for both methods, it would have been possible to increase the 
number of configurations or the number of test particles in Widoms method by a factor of 
10. Increasing the number of test particles would have lead to no improvements, as their 
number is already very high. Increasing the number of configurations by a factor of 10 
would have lead to a reduction of the statistical uncertainties by about a factor of 3 at 
best. 

For the 2CLJQ fluid (L* = 0.8, Q*^ = 4) at the lower temperature {T/T^ ^ 0.55) a study 
of the influence of the free parameters of the gradual insertion method on the chemical 
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potential is given in Table 2. The reference point is Mq—IO, Mp—lO, and Mp— 50. These 
parameters were altered up and down by a factor of 2.5. The increase of the number of 
attempts to change the state of the fluctuating particle from Mc=4 to Mc— 10 yields 
lower uncertainties, whereas upon further increase to Mc=25 no benefit was observed. 
The variation of the number of attempts to translate or rotate the fluctuating particle 
(Mp) seems to have no influence on the uncertainty in the investigated range of that 
parameter. The number of attempts to find a host particle for preferential translation 
or rotation (Mp) has the clearest influence on the uncertainty. Note that an increase 
above MpfalOO is not recommended. For instance, at Mp— 125 we have a relation of 125 
preferential translation moves to one ordinary translation move. This extreme ratio leads 
to the slight deviation shown in the last line of Table 2. 

5 Conclusion 

This investigation shows that NpT simulations with gradual insertion combined with 
preferential sampling can yield results for the chemical potential of realistic strongly in- 
teracting molecular fluids at high densities with clearly improved statistics. It can be 
applied at state points where conventional test particle insertion breaks down. 
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Table 1 

Density, configurational energy, and chemical potential for different 2CLJQ fluids close to their 
bubble lines at T/Tc ~ 0.55 (upper blocks) and T/Tc ~ 0.8 (lower blocks). Method A: gradual 
insertion with the parameters Mc = 10, Mp = 10, and Mp = 50; Method B: Widom's test 
particle method [21]. The numbers in paranthesis denote the uncertainty in the last digits. 





p* 


P* 


u* 


fi/kT 


Method 


L* = 0.4, Q*2 = 


1.7875 





0.5960 (3) 


-15.524 (8) 


-6.31 (2) 


A 


0.5952 (2) 


-15.504 (6) 


-6.34 (6) 


B 


2.6 


0.06 


0.4747 (8) 


-11.80 (2) 


-3.867 (4) 


A 


0.4718 (6) 


-11.74 (1) 


-3.882 (7) 


B 


L* = 0.4, Q*2 = 4 


2.09 





0.6566 (3) 


-23.58 (1) 


-7.05 (3) 


A 


0.6562 (2) 


-23.57 (1) 


-6.8 (5) 


B 


3.04 


0.07 


0.5180 (8) 


-16.70 (3) 


-3.95 (1) 


A 


0.5154 (4) 


-16.61 (2) 


-3.98 (1) 


B 


L* = 0.6, Q*2 = 


1.4025 





0.5008 (2) 


-12.706 (6) 


-6.59 (1) 


A 


0.5001 (2) 


-12.689 (5) 


-6.57 (7) 


B 


2.04 


0.04 


0.3947 (7) 


-9.51 (2) 


-4.03 (1) 


A 


0.3931 (4) 


-9.47 (1) 


-4.058 (8) 


B 


L* = 0.6, Q*2 = 4 


1.628 





0.5521 (2) 


-19.06 (1) 


-7.31 (3) 


A 


0.5513 (2) 


-19.03 (1) 


-7.6 (4) 


B 


2.368 


0.04 


0.4310 (7) 


-13.33 (2) 


-4.16 (1) 


A 


0.4298 (5) 


-13.29 (2) 


-4.16 (1) 


B 
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Table 1 - continued 





P* 


P* 


u* 


li/kT 


Method 


L* = 0.8, Q*2 = 


1.177 





0.4389 (2) 


-11.124 (4) 


-6.82 (1) 


A 


0.4384 (1) 


-11.109 (4) 


-6.7 (1) 


B 


1.712 


0.03 


0.3421 (6) 


-8.19 (1) 


-4.161 (8) 


A 


0.3398 (5) 


-8.14 (1) 


-4.192 (8) 


B 


L* = 0.8, Q*2 = 4 


1.342 





0.49G3 (2) 


-17.31 (1) 


-7.91 (2) 


A 


0.4955 (2) 


-17.276 (8) 


-8.3 (8) 


B 


1.952 


0.03 


0.3932 (4) 


-12.36 (1) 


-4.495 (8) 


A 


0.3928 (4) 


-12.35 (1) 


-4.51 (2) 


B 
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Table 2 

Influence of the gradual insertion parameters on the chemical potential of the 2CLJQ fluid 
(L* = 0.8, Q*"^ = 4) at T* = 1.342, p* = 0. The numbers in paranthesis denote the uncertainty 
in the last digits. 



Mc 


Mf 


Mp 


fi/kT 


10 


10 


50 


-7.97 (2) 


4 


10 


50 


-8.00 (5) 


25 


10 


50 


-8.04 (3) 


10 


4 


50 


-8.00 (1) 


10 


25 


50 


-7.99 (4) 


10 


10 


20 


-7.97 (6) 


10 


10 


125 


-7.94 (1) 
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List of Figures 



1 Running averages of tlie cliemical potential over Monte-Carlo loops, or 
respectively, molecular dynamics time steps of independent simulation 
runs for the 2CLJQ fluid {L* = 0.6, Q*^ = 4) at T* = 1.628, p* = using 
N = 512 particles and a cut-off radius = 4 o". Thick line: Monte-Carlo 
simulation using gradual insertion with the parameters Mc=10, Mf=10, 
and Mp=50. Thin lines: Three different molecular dynamics simulation 
runs using Widom's method with 2 000 test particles at each time step. 
Dashed lines: straight horizontal guide to the eye. 
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Fig. 1. 
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